Bayesian GLM Part4

Preparations

Load the necessary libraries

library(car)
library(rstanarm) # for fitting models in STAN
library(brms) # for fitting models in STAN
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(DHARMa) # for residual diagnostics
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(ggeffects) # for partial effects plots
library(tidyverse) # for data wrangling etc
library(patchwork) # for multiple plots

Scenario

Loyn (1987) modelled the abundance of forest birds with six predictor variables (patch area, distance to nearest patch, distance to nearest larger patch, grazing intensity, altitude and years since the patch had been isolated).

Regent honeyeater

Format of loyn.csv data file

ABUND	DIST	LDIST	AREA	GRAZE	ALT	YR.ISOL
..	..	..	..	..	..	..

ABUND	Abundance of forest birds in patch- response variable
DIST	Distance to nearest patch - predictor variable
LDIST	Distance to nearest larger patch - predictor variable
AREA	Size of the patch - predictor variable
GRAZE	Grazing intensity (1 to 5, representing light to heavy) - predictor variable
ALT	Altitude - predictor variable
YR.ISOL	Number of years since the patch was isolated - predictor variable

The aim of the analysis is to investigate the effects of a range of predictors on the abundance of forest birds.

Read in the data

loyn <- read_csv("../public/data/loyn.csv", trim_ws = TRUE)
glimpse(loyn)

## Rows: 56
## Columns: 7
## $ ABUND   <dbl> 5.3, 2.0, 1.5, 17.1, 13.8, 14.1, 3.8, 2.2, 3.3, 3.0, 27.6, 1.8…
## $ AREA    <dbl> 0.1, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.…
## $ YR.ISOL <dbl> 1968, 1920, 1900, 1966, 1918, 1965, 1955, 1920, 1965, 1900, 19…
## $ DIST    <dbl> 39, 234, 104, 66, 246, 234, 467, 284, 156, 311, 66, 93, 39, 40…
## $ LDIST   <dbl> 39, 234, 311, 66, 246, 285, 467, 1829, 156, 571, 332, 93, 39, …
## $ GRAZE   <dbl> 2, 5, 5, 3, 5, 3, 5, 5, 4, 5, 3, 5, 2, 1, 5, 5, 3, 3, 3, 2, 2,…
## $ ALT     <dbl> 160, 60, 140, 160, 140, 130, 90, 60, 130, 130, 210, 160, 210, …

Exploratory data analysis

When we explored this analysis from a frequentist perspective, we decided on a log-normal like model. This was a model that was fit against a Gaussian distribution, yet with a log-link. We will replicate that model here in a Bayesian framework.

In the previous exploration of this model, we elected to treat Grazing intensity as a categorical variable - we will again code Grazing intensity as a categorical variable.

loyn <- loyn %>% mutate(fGRAZE = factor(GRAZE))

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ log(\mu_i) = \boldsymbol{\beta} \bf{X_i} \]

where \(\boldsymbol{\beta}\) is a vector of effects parameters and \(\bf{X}\) is a model matrix representing the additive effects of the scaled versions of distance (ln), distance to the nearest large patch (ln), patch area (ln), grazing intensity, year of isolation and altitude on the abundance of forest birds.

To re-acquaint ourselves with the data, I we will revisit the scatterplot matrix that we generated prior to the frequentist analysis.

scatterplotMatrix(~ ABUND + DIST + LDIST + AREA + GRAZE + ALT + YR.ISOL,
  data = loyn,
  diagonal = list(method = "boxplot")
)

We again notice that DIST, LDIST and AREA are skewed, so we will normalise them via a logarithmic transformation.

scatterplotMatrix(~ ABUND + log(DIST) + log(LDIST) + log(AREA) + GRAZE + ALT + YR.ISOL,
  data = loyn,
  diagonal = list(method = "boxplot")
)

Fit the model

frequentist

loyn.glm <- glm(ABUND ~ scale(log(DIST)) + scale(log(LDIST)) + scale(log(AREA)) +
  fGRAZE + scale(ALT) + scale(YR.ISOL),
data = loyn,
family = gaussian(link = "log")
)
loyn.glm %>% summary()

## 
## Call:
## glm(formula = ABUND ~ scale(log(DIST)) + scale(log(LDIST)) + 
##     scale(log(AREA)) + fGRAZE + scale(ALT) + scale(YR.ISOL), 
##     family = gaussian(link = "log"), data = loyn)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -15.9887   -3.3663   -0.7579    4.0166   12.8221  
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.092248   0.112856  27.400  < 2e-16 ***
## scale(log(DIST))  -0.018067   0.057247  -0.316  0.75373    
## scale(log(LDIST))  0.057086   0.059984   0.952  0.34623    
## scale(log(AREA))   0.203976   0.059620   3.421  0.00132 ** 
## fGRAZE2            0.039644   0.148978   0.266  0.79134    
## fGRAZE3            0.019654   0.137752   0.143  0.88717    
## fGRAZE4           -0.001197   0.156199  -0.008  0.99392    
## fGRAZE5           -1.121563   0.343631  -3.264  0.00208 ** 
## scale(ALT)        -0.003237   0.048607  -0.067  0.94719    
## scale(YR.ISOL)    -0.018246   0.074404  -0.245  0.80737    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 42.40928)
## 
##     Null deviance: 6337.9  on 55  degrees of freedom
## Residual deviance: 1950.8  on 46  degrees of freedom
## AIC: 379.76
## 
## Number of Fisher Scoring iterations: 6

rstanarm

Using default priors

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.

loyn.rstanarm <- stan_glm(ABUND ~ scale(log(DIST), scale = FALSE) +
  scale(log(LDIST), scale = FALSE) +
  scale(log(AREA), scale = FALSE) +
  fGRAZE +
  scale(ALT, scale = FALSE) +
  scale(YR.ISOL, scale = FALSE),
data = loyn,
family = gaussian(link = "log"),
iter = 5000, warmup = 1000,
chains = 3, thin = 5, refresh = 0
)

Conclusions:

• the warning messages suggest that this model has not performed well
  • there are a large number of divergent transitions. Divergent transitions occur when the sampler is in a ‘sharp’ part of the posterior and the step length is too large causing the sampler to shoot off the posterior and reject the sample. This leads to inefficient and ineffective MCMC sampling. To alleviate this, we can either:
    • review the model itself - it might be misspecified
    • adjust the adaptive delta - a parameter that governs the degree of step-length learning that occurs during the warmup stage - the more it learns the less likely the sampler will be to overshoot (although it will take longer to learn)
    • review the priors
  • the largest R-hat value is large. This suggests that the chains have not converged well
  • the effective sample sizes is too low, suggesting that there might be issues with the priors.

prior_summary(loyn.rstanarm)

## Priors for model 'loyn.rstanarm' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 0, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0, scale = 27)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [28.17,20.27,14.35,...])
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.093)
## ------
## See help('prior_summary.stanreg') for more details

Conclusions:

• the default priors appear to be overly wide. It seems that the adjusted priors are neglecting the log link here.

Assessing priors

loyn.rstanarm1 <- update(loyn.rstanarm, prior_PD = TRUE)

ggemmeans(loyn.rstanarm1, ~AREA) %>% plot(add.data = TRUE) + scale_y_log10()

Conclusions:

• we see that the range of predictions are very wide and the slope could range from strongly negative to strongly positive.

Defining priors

• \(\beta_0\): normal centred at 3 with a standard deviation of 10 (on the log scale, these are the equivalent of 20.1 and 22026 respectively)
• \(\beta_1\): normal centred at 0 with a standard deviation of 2.5 (again, consider what this would be on a response scale)
• \(\sigma\): half cauchy (flat Gaussian) prior with a scale of 5.

I will also overlay the raw data for comparison.

loyn.rstanarm2 <- stan_glm(ABUND ~ scale(log(DIST)) +
  scale(log(LDIST)) +
  scale(log(AREA)) +
  fGRAZE +
  scale(ALT) +
  scale(YR.ISOL),
data = loyn,
family = gaussian(link = "log"),
prior_intercept = normal(3, 10, autoscale = FALSE),
prior = normal(0, 2.5, autoscale = FALSE),
prior_aux = cauchy(0, 5),
prior_PD = TRUE,
iter = 5000, thin = 5, chains = 3, warmup = 2000,
refresh = 0
)

Conclusions:

• there are no longer any warning messages.

ggemmeans(loyn.rstanarm2, ~AREA) %>%
  plot(add.data = TRUE) + scale_y_log10()

Now lets refit, conditioning on the data.

loyn.rstanarm3 <- update(loyn.rstanarm2, prior_PD = FALSE)

Plotting prior and posterior

posterior_vs_prior(loyn.rstanarm3,
  color_by = "vs", group_by = TRUE,
  facet_args = list(scales = "free_y")
)

Conclusions:

• in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

# ggemmeans(loyn.rstanarm3,  ~AREA) %>% plot(add.data=TRUE) + scale_y_log10()

brms

Using default priors

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

loyn.brm <- brm(bf(ABUND ~ scale(log(DIST)) +
  scale(log(LDIST)) +
  scale(log(AREA)) +
  fGRAZE +
  scale(ALT) +
  scale(YR.ISOL),
family = lognormal()
),
data = loyn,
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)

loyn.brm %>% prior_summary()

##                 prior     class          coef group resp dpar nlpar bound       source
##                (flat)         b                                                default
##                (flat)         b       fGRAZE2                             (vectorized)
##                (flat)         b       fGRAZE3                             (vectorized)
##                (flat)         b       fGRAZE4                             (vectorized)
##                (flat)         b       fGRAZE5                             (vectorized)
##                (flat)         b      scaleALT                             (vectorized)
##                (flat)         b  scalelogAREA                             (vectorized)
##                (flat)         b  scalelogDIST                             (vectorized)
##                (flat)         b scalelogLDIST                             (vectorized)
##                (flat)         b  scaleYR.ISOL                             (vectorized)
##  student_t(3, 3, 2.5) Intercept                                                default
##  student_t(3, 0, 2.5)     sigma                                                default

Assessing priors

priors <- prior(normal(0, 2.5), class = "b")

loyn.brm1 <- brm(bf(ABUND ~ scale(log(DIST)) +
  scale(log(LDIST)) +
  scale(log(AREA)) +
  fGRAZE +
  scale(ALT) +
  scale(YR.ISOL),
family = gaussian(link = "log")
),
data = loyn,
prior = priors,
sample_prior = "only",
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)

## Individual plots
loyn.brm1 %>%
  ggemmeans(~AREA) %>%
  plot(add.data = TRUE) + scale_y_log10()

## All effects
loyn.brm1 %>%
  conditional_effects() %>%
  plot(points = TRUE, ask = FALSE, plot = FALSE) %>%
  wrap_plots()

## All effects log y axis
## Do above, but then modify each list item
loyn.brm1 %>%
  conditional_effects() %>%
  plot(points = TRUE, ask = FALSE, plot = FALSE) %>%
  lapply(function(x) x + scale_y_log10()) %>%
  wrap_plots()

Conclusions:

• we see that the range of predictions is fairly wide and the slope could range from strongly negative to strongly positive.

Defining priors

priors <- prior(normal(0, 10), class = "Intercept") +
  prior(normal(0, 1), class = "b") +
  prior(gamma(2, 1), class = "sigma")

loyn.brm2 <- brm(bf(ABUND ~ scale(log(DIST)) +
  scale(log(LDIST)) +
  scale(log(AREA)) +
  fGRAZE +
  scale(ALT) +
  scale(YR.ISOL),
family = gaussian(link = "log")
),
data = loyn,
prior = priors,
sample_prior = "only",
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)

loyn.brm2 %>%
  ggemmeans(~AREA) %>%
  plot(add.data = TRUE) +
  scale_y_log10()

loyn.brm2 %>%
  conditional_effects() %>%
  plot(points = TRUE, ask = FALSE, plot = FALSE) %>%
  lapply(function(x) x + scale_y_log10()) %>%
  wrap_plots()

loyn.brm3 <- update(loyn.brm2, sample_prior = "yes", refresh = 0)

Plotting prior and posterior

loyn.brm3 %>% get_variables()

##  [1] "b_Intercept"     "b_scalelogDIST"  "b_scalelogLDIST" "b_scalelogAREA" 
##  [5] "b_fGRAZE2"       "b_fGRAZE3"       "b_fGRAZE4"       "b_fGRAZE5"      
##  [9] "b_scaleALT"      "b_scaleYR.ISOL"  "sigma"           "prior_Intercept"
## [13] "prior_b"         "prior_sigma"     "lp__"            "accept_stat__"  
## [17] "stepsize__"      "treedepth__"     "n_leapfrog__"    "divergent__"    
## [21] "energy__"

loyn.brm3 %>%
  hypothesis("Intercept = 0", class = "b") %>%
  plot()

loyn.brm3 %>%
  hypothesis("Intercept = 0", class = "prior") %>%
  plot()

loyn.brm3 %>%
  hypothesis("scalelogAREA = 0") %>%
  plot()

loyn.brm3 %>%
  hypothesis("sigma = 0", class = "") %>%
  plot()

loyn.brm3 %>%
  posterior_samples() %>%
  select(-`lp__`) %>%
  pivot_longer(everything(), names_to = "key") %>%
  mutate(
    Type = ifelse(str_detect(key, "prior"), "Prior", "b"),
    Class = ifelse(str_detect(key, "Intercept"), "Intercept",
      ifelse(str_detect(key, "b"), "b", "sigma")
    ),
    Par = str_replace(key, "b_", "")
  ) %>%
  ggplot(aes(x = Type, y = value, color = Par)) +
  stat_pointinterval(position = position_dodge()) +
  facet_wrap(~Class, scales = "free")

Exploring the stan code

loyn.brm3 %>% standata()

## $N
## [1] 56
## 
## $Y
##  [1]  5.3  2.0  1.5 17.1 13.8 14.1  3.8  2.2  3.3  3.0 27.6  1.8 21.2 14.6  8.0
## [16]  3.5 29.0  2.9 24.3 19.4 24.4  5.0 15.8 25.3 19.5 20.9 16.3 18.8 19.9 13.0
## [31]  6.8 21.7 27.8 26.8 16.6 30.4 11.5 26.0 25.7 12.7 23.5 24.9 29.0 28.3 28.3
## [46] 32.0 37.7 39.6 29.6 31.0 34.4 27.3 30.5 33.0 29.5 30.9
## 
## $K
## [1] 10
## 
## $X
##    Intercept scalelogDIST scalelogLDIST scalelogAREA fGRAZE2 fGRAZE3 fGRAZE4
## 1          1  -1.51019511   -1.65892116  -2.37802367       1       0       0
## 2          1   0.37029448   -0.30532977  -1.51765965       0       0       0
## 3          1  -0.48079400   -0.09042449  -1.51765965       0       0       0
## 4          1  -0.95804924   -1.26148217  -1.14712103       0       1       0
## 5          1   0.42278148   -0.26754922  -1.14712103       0       0       0
## 6          1   0.37029448   -0.15637842  -1.14712103       0       1       0
## 7          1   1.09552219    0.21669491  -1.14712103       0       0       0
## 8          1   0.57353754    1.24803687  -1.14712103       0       0       0
## 9          1  -0.05524976   -0.61163991  -1.14712103       0       0       1
## 10         1   0.66885367    0.36858640  -1.14712103       0       0       0
## 11         1  -0.95804924   -0.04106159  -0.77658241       0       1       0
## 12         1  -0.59812145   -1.00240327  -0.77658241       0       0       0
## 13         1  -1.51019511   -1.65892116  -0.77658241       1       0       0
## 14         1   0.93822292    0.10346964  -0.77658241       0       0       0
## 15         1   0.47682817   -0.22864597  -0.77658241       0       0       0
## 16         1  -0.24660010    0.43442972  -0.77658241       0       0       0
## 17         1  -1.93573935   -1.96523129  -0.55983122       0       1       0
## 18         1  -1.93573935   -1.96523129  -0.55983122       0       1       0
## 19         1  -1.48362354   -1.63979473  -0.40604380       0       1       0
## 20         1   0.47682817   -0.22864597  -0.40604380       1       0       0
## 21         1   0.37029448    0.29645165  -0.40604380       1       0       0
## 22         1  -1.93573935    1.38927509  -0.40604380       0       0       0
## 23         1  -1.51019511   -1.65892116  -0.28675700       0       1       0
## 24         1   0.85682391    0.04487798  -0.28675700       0       0       0
## 25         1  -0.59812145   -0.33160908  -0.18929260       0       1       0
## 26         1  -0.03525827    0.79868571  -0.18929260       0       1       0
## 27         1   0.57722655    0.69705984  -0.10688762       0       1       0
## 28         1  -0.22265561   -0.73213997  -0.10688762       0       0       1
## 29         1   0.50481706   -0.20849935  -0.03550518       0       0       1
## 30         1   0.79284436    1.26397616   0.02745860       0       0       0
## 31         1   0.75311975    0.29645165   0.08378161       0       1       0
## 32         1  -1.89613008   -1.93672022   0.13473198       0       0       1
## 33         1  -0.03525827   -0.59724988   0.18124602       0       0       1
## 34         1  -0.22265561   -0.73213997   0.18124602       0       0       1
## 35         1  -0.12477989   -0.66168825   0.22403479       1       0       0
## 36         1  -0.12477989    0.09591504   0.30053281       0       1       0
## 37         1   0.90372228    1.51230754   0.36744180       0       0       0
## 38         1  -1.48362354    1.66778539   0.39799721       1       0       0
## 39         1   0.50481706    0.99128245   0.42690016       0       0       1
## 40         1   0.66885367    1.53439756   0.50527059       0       0       0
## 41         1   1.35327186    0.40222517   0.59457340       0       0       0
## 42         1   1.25762760    0.59446805   0.65294853       0       1       0
## 43         1   0.24667868   -0.39430941   0.70557205       0       0       0
## 44         1  -0.95804924   -0.01204503   0.73798041       0       0       0
## 45         1   0.57722655   -0.51197475   0.82485885       1       0       0
## 46         1  -0.59812145   -1.00240327   0.87580921       0       1       0
## 47         1   0.47682817    0.98837574   0.92232325       0       1       0
## 48         1   2.26783778    1.12640241   0.93334579       0       0       0
## 49         1   0.92772762    1.07832552   0.94414564       0       0       0
## 50         1   1.09552219    0.21669491   1.01418997       0       0       0
## 51         1  -1.51019511    0.29645165   1.29286187       1       0       0
## 52         1   0.93822292    1.91565164   1.35582564       0       0       0
## 53         1   1.09552219    1.39201101   1.47113789       0       0       0
## 54         1  -0.12477989   -0.07123733   1.50961306       0       0       0
## 55         1   0.75311975    1.00336997   2.53095496       0       0       0
## 56         1   0.73743154   -0.04106159   2.85111978       0       0       0
##    fGRAZE5    scaleALT scaleYR.ISOL
## 1        0  0.31588146    0.7134084
## 2        1 -1.98143824   -1.1629534
## 3        1 -0.14358248   -1.9447708
## 4        0  0.31588146    0.6352266
## 5        1 -0.14358248   -1.2411351
## 6        0 -0.37331445    0.5961358
## 7        1 -1.29224233    0.2052271
## 8        1 -1.98143824   -1.1629534
## 9        0 -0.37331445    0.5961358
## 10       1 -0.37331445   -1.9447708
## 11       0  1.46454131   -0.9284082
## 12       1  0.31588146   -2.3356795
## 13       0  1.46454131    0.9088627
## 14       0  1.46454131    0.8697719
## 15       1 -0.60304642   -1.9447708
## 16       1 -0.02871650   -1.9447708
## 17       0 -0.83277839    0.4788632
## 18       0 -0.14358248    0.5961358
## 19       0  1.00507737    0.4006814
## 20       0 -1.29224233    0.1270453
## 21       0  1.69427328    0.9088627
## 22       1 -0.60304642   -1.0456808
## 23       0 -0.37331445    0.5961358
## 24       0 -1.06251036    0.6743175
## 25       0  0.54561343   -2.3356795
## 26       0  0.08614949    0.4006814
## 27       0 -0.37331445    0.5961358
## 28       0  1.46454131    0.4006814
## 29       0 -0.60304642    0.9088627
## 30       1 -1.29224233   -1.5538621
## 31       0 -0.83277839    0.4788632
## 32       0  0.66047941    0.4006814
## 33       0 -0.83277839    0.5179540
## 34       0  1.23480934    0.4006814
## 35       0  1.00507737    0.7134084
## 36       0 -0.60304642    0.6352266
## 37       1 -1.06251036   -1.1629534
## 38       0  1.00507737    0.6352266
## 39       0  0.08614949    0.9088627
## 40       1 -1.29224233   -1.2411351
## 41       0 -0.14358248    0.5179540
## 42       0 -0.37331445    0.5961358
## 43       0  1.00507737    0.9479536
## 44       0 -0.83277839    0.5961358
## 45       0 -0.60304642    0.4788632
## 46       0  1.00507737    0.4006814
## 47       0 -0.60304642   -0.8502264
## 48       0  0.77534540    0.8697719
## 49       0 -0.14358248    0.6743175
## 50       0  0.43074744    0.5179540
## 51       0  0.66047941    1.0261353
## 52       0 -1.75170627    0.5570449
## 53       0  0.31588146    0.5570449
## 54       0  1.00507737   -0.3811360
## 55       0  1.00507737    0.7915901
## 56       0  2.61320116   -0.6547721
## attr(,"assign")
##  [1] 0 1 2 3 4 4 4 4 5 6
## attr(,"contrasts")
## attr(,"contrasts")$fGRAZE
##   2 3 4 5
## 1 0 0 0 0
## 2 1 0 0 0
## 3 0 1 0 0
## 4 0 0 1 0
## 5 0 0 0 1
## 
## 
## $prior_only
## [1] 0
## 
## attr(,"class")
## [1] "standata"

loyn.brm3 %>% stancode()

## // generated with brms 2.14.4
## functions {
## }
## data {
##   int<lower=1> N;  // total number of observations
##   vector[N] Y;  // response variable
##   int<lower=1> K;  // number of population-level effects
##   matrix[N, K] X;  // population-level design matrix
##   int prior_only;  // should the likelihood be ignored?
## }
## transformed data {
##   int Kc = K - 1;
##   matrix[N, Kc] Xc;  // centered version of X without an intercept
##   vector[Kc] means_X;  // column means of X before centering
##   for (i in 2:K) {
##     means_X[i - 1] = mean(X[, i]);
##     Xc[, i - 1] = X[, i] - means_X[i - 1];
##   }
## }
## parameters {
##   vector[Kc] b;  // population-level effects
##   real Intercept;  // temporary intercept for centered predictors
##   real<lower=0> sigma;  // residual SD
## }
## transformed parameters {
## }
## model {
##   // likelihood including all constants
##   if (!prior_only) {
##     // initialize linear predictor term
##     vector[N] mu = Intercept + Xc * b;
##     for (n in 1:N) {
##       // apply the inverse link function
##       mu[n] = exp(mu[n]);
##     }
##     target += normal_lpdf(Y | mu, sigma);
##   }
##   // priors including all constants
##   target += normal_lpdf(b | 0, 1);
##   target += normal_lpdf(Intercept | 0, 10);
##   target += gamma_lpdf(sigma | 2, 1);
## }
## generated quantities {
##   // actual population-level intercept
##   real b_Intercept = Intercept - dot_product(means_X, b);
##   // additionally draw samples from priors
##   real prior_b = normal_rng(0,1);
##   real prior_Intercept = normal_rng(0,10);
##   real prior_sigma = gamma_rng(2,1);
##   // use rejection sampling for truncated priors
##   while (prior_sigma < 0) {
##     prior_sigma = gamma_rng(2,1);
##   }
## }

MCMC sampling diagnostics

In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.

There are a wide variety of tests that range from the big picture, overall chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.

rstanarm

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()

## bayesplot MCMC module:
##   mcmc_acf
##   mcmc_acf_bar
##   mcmc_areas
##   mcmc_areas_data
##   mcmc_areas_ridges
##   mcmc_areas_ridges_data
##   mcmc_combo
##   mcmc_dens
##   mcmc_dens_chains
##   mcmc_dens_chains_data
##   mcmc_dens_overlay
##   mcmc_hex
##   mcmc_hist
##   mcmc_hist_by_chain
##   mcmc_intervals
##   mcmc_intervals_data
##   mcmc_neff
##   mcmc_neff_data
##   mcmc_neff_hist
##   mcmc_nuts_acceptance
##   mcmc_nuts_divergence
##   mcmc_nuts_energy
##   mcmc_nuts_stepsize
##   mcmc_nuts_treedepth
##   mcmc_pairs
##   mcmc_parcoord
##   mcmc_parcoord_data
##   mcmc_rank_hist
##   mcmc_rank_overlay
##   mcmc_recover_hist
##   mcmc_recover_intervals
##   mcmc_recover_scatter
##   mcmc_rhat
##   mcmc_rhat_data
##   mcmc_rhat_hist
##   mcmc_scatter
##   mcmc_trace
##   mcmc_trace_data
##   mcmc_trace_highlight
##   mcmc_violin

Of these, we will focus on:

• mcmc_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

plot(loyn.rstanarm3, plotfun = "mcmc_trace")

The chains appear well mixed and very similar

• acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

plot(loyn.rstanarm3, "acf_bar")

There is no evidence of auto-correlation in the MCMC samples

• Rhat: Rhat is a measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

plot(loyn.rstanarm3, "rhat_hist")

All Rhat values are below 1.05, suggesting the chains have converged.

• neff (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

plot(loyn.rstanarm3, "neff_hist")

Ratios all very high.

More diagnostics

plot(loyn.rstanarm3, "combo")

plot(loyn.rstanarm3, "violin")

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

stan_trace(loyn.rstanarm3)

The chains appear well mixed and very similar

• stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

stan_ac(loyn.rstanarm3)

There is no evidence of auto-correlation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

stan_rhat(loyn.rstanarm3)

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

stan_ess(loyn.rstanarm3)

Ratios all very high.

stan_dens(loyn.rstanarm3, separate_chains = TRUE)

ggmcmc

The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

loyn.ggs <- ggs(loyn.rstanarm3)
ggs_traceplot(loyn.ggs)

The chains appear well mixed and very similar

• gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

ggs_autocorrelation(loyn.ggs)

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

ggs_Rhat(loyn.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

ggs_effective(loyn.ggs)

Ratios all very high.

More diagnostics

ggs_crosscorrelation(loyn.ggs)

ggs_grb(loyn.ggs)

brm

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()

## bayesplot MCMC module:
##   mcmc_acf
##   mcmc_acf_bar
##   mcmc_areas
##   mcmc_areas_data
##   mcmc_areas_ridges
##   mcmc_areas_ridges_data
##   mcmc_combo
##   mcmc_dens
##   mcmc_dens_chains
##   mcmc_dens_chains_data
##   mcmc_dens_overlay
##   mcmc_hex
##   mcmc_hist
##   mcmc_hist_by_chain
##   mcmc_intervals
##   mcmc_intervals_data
##   mcmc_neff
##   mcmc_neff_data
##   mcmc_neff_hist
##   mcmc_nuts_acceptance
##   mcmc_nuts_divergence
##   mcmc_nuts_energy
##   mcmc_nuts_stepsize
##   mcmc_nuts_treedepth
##   mcmc_pairs
##   mcmc_parcoord
##   mcmc_parcoord_data
##   mcmc_rank_hist
##   mcmc_rank_overlay
##   mcmc_recover_hist
##   mcmc_recover_intervals
##   mcmc_recover_scatter
##   mcmc_rhat
##   mcmc_rhat_data
##   mcmc_rhat_hist
##   mcmc_scatter
##   mcmc_trace
##   mcmc_trace_data
##   mcmc_trace_highlight
##   mcmc_violin

Of these, we will focus on:

• mcmc_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

loyn.brm3 %>% mcmc_plot(type = "trace")

The chains appear well mixed and very similar

• acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

loyn.brm3 %>% mcmc_plot(type = "acf_bar")

There is no evidence of autocorrelation in the MCMC samples

• Rhat: Rhat is a measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

loyn.brm3 %>% mcmc_plot(type = "rhat_hist")

All Rhat values are below 1.05, suggesting the chains have converged.

• neff (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

loyn.brm3 %>% mcmc_plot(type = "neff_hist")

Ratios all very high.

More diagnostics

loyn.brm3 %>% mcmc_plot(type = "combo")

loyn.brm3 %>% mcmc_plot(type = "violin")

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

loyn.brm3$fit %>% stan_trace()

The chains appear well mixed and very similar

• stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

loyn.brm3$fit %>% stan_ac()

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

loyn.brm3$fit %>% stan_rhat()

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

loyn.brm3$fit %>% stan_ess()

Ratios all very high.

loyn.brm3$fit %>% stan_dens(separate_chains = TRUE)

ggmcmc

The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

loyn.ggs <- loyn.brm3 %>% ggs(inc_warmup = FALSE, burnin = FALSE)
loyn.ggs %>% ggs_traceplot()

The chains appear well mixed and very similar

• gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

loyn.ggs %>% ggs_autocorrelation()

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

loyn.ggs %>% ggs_Rhat()

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

loyn.ggs %>% ggs_effective()

Ratios all very high.

More diagnostics

loyn.ggs %>% ggs_crosscorrelation()

loyn.ggs %>% ggs_grb()

Model validation

rstanarm

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_dens_overlay_grouped
##   ppc_ecdf_overlay
##   ppc_ecdf_overlay_grouped
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_km_overlay
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_data
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

• dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

pp_check(loyn.rstanarm3, plotfun = "dens_overlay")

The model draws appear deviate from the observed data.

• error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.

pp_check(loyn.rstanarm3, plotfun = "error_scatter_avg")

The predictive error seems to be related to the predictor - the model performs poorest at higher mussel clump areas.

• error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such. Again, this is not interpretable for binary data.

pp_check(loyn.rstanarm3, x = loyn$AREA, plotfun = "error_scatter_avg_vs_x")

• intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

pp_check(loyn.rstanarm3, x = loyn$AREA, plotfun = "intervals")

The modelled predictions seem to underestimate the uncertainty with increasing mussel clump area.

• ribbon: this is just an alternative way of expressing the above plot.

pp_check(loyn.rstanarm3, x = loyn$AREA, plotfun = "ribbon")

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(loyn.rstanarm3)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

• simulated (predicted) responses associated with each observation.
• observed values
• fitted (predicted) responses (averaged) associated with each observation

preds <- posterior_predict(loyn.rstanarm3, nsamples = 250, summary = FALSE)
loyn.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = loyn$ABUND,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = FALSE
)
plot(loyn.resids)

Conclusions:

• the simulated residuals suggest a general lack of fit due to over dispersion and outliers
• perhaps we should explore a negative binomial model

brms

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_dens_overlay_grouped
##   ppc_ecdf_overlay
##   ppc_ecdf_overlay_grouped
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_km_overlay
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_data
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

• dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

loyn.brm3 %>% pp_check(type = "dens_overlay")

The model draws appear deviate from the observed data.

• error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.

loyn.brm3 %>% pp_check(type = "error_scatter_avg")

The predictive error seems to be related to the predictor - the model performs poorest at higher mussel clump areas.

• error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such. Again, this is not interpretable for binary data.

loyn.brm3 %>% pp_check(x = "AREA", type = "error_scatter_avg_vs_x")

• intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

loyn.brm3 %>% pp_check(x = "AREA", type = "intervals")

The modelled predictions seem to underestimate the uncertainty with increasing mussel clump area.

• ribbon: this is just an alternative way of expressing the above plot.

loyn.brm3 %>% pp_check(x = "AREA", type = "ribbon")

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(loyn.brm3)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

• simulated (predicted) responses associated with each observation.
• observed values
• fitted (predicted) responses (averaged) associated with each observation

preds <- loyn.brm3 %>% posterior_predict(nsamples = 250, summary = FALSE)
loyn.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = loyn$ABUND,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = FALSE
)
loyn.resids %>% plot()

Conclusions:

• the simulated residuals suggest a general lack of fit due to over-dispersion and outliers
• perhaps we should explore a negative binomial model

Partial effects plots

rstanarm

ggpredict

loyn.rstanarm3 %>%
  ggpredict() %>%
  plot(add.data = TRUE, facet = TRUE)

ggemmeans

# loyn.rstanarm3 %>% ggemmeans(~AREA,  type='fixed') %>% plot(add.data=TRUE) + scale_y_log10()

fitted_draws

loyn.rstanarm3 %>%
  fitted_draws(newdata = loyn) %>%
  median_hdci() %>%
  ggplot(aes(x = AREA, y = .value)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper), fill = "blue", alpha = 0.3) +
  geom_line() +
  geom_point(data = loyn, aes(y = ABUND, x = AREA)) +
  scale_y_log10()

brms

conditional_effects

loyn.brm3 %>%
  conditional_effects() %>%
  plot(ask = FALSE, points = TRUE, plot = FALSE) %>%
  wrap_plots()

g <- loyn.brm3 %>%
  conditional_effects() %>%
  plot(ask = FALSE, points = TRUE, plot = FALSE)
library(patchwork)
length(g)

## [1] 6

(g[[1]] + scale_x_log10()) +
  (g[[2]] + scale_x_log10()) +
  (g[[3]] + scale_x_log10()) +
  g[[4]] +
  g[[5]] +
  g[[6]]

ggpredict

loyn.brm3 %>%
  ggpredict() %>%
  plot(add.data = TRUE) %>%
  wrap_plots()

ggemmeans

## loyn.brm3 %>%
##     ggemmeans(~AREA) %>%
##     plot(add.data=TRUE) %>%
##     wrap_plots()

fitted_draws

loyn.brm3 %>%
  fitted_draws(newdata = loyn) %>%
  median_hdci() %>%
  ggplot(aes(x = AREA, y = .value, colour = fGRAZE, fill = fGRAZE)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper), colour = NA, alpha = 0.3) +
  geom_line() +
  geom_point(data = loyn, aes(y = ABUND, x = AREA)) +
  scale_y_log10() +
  scale_x_log10()

Model investigation

rstanarm

rstanarm captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

summary(loyn.rstanarm3)

## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [log]
##  formula:      ABUND ~ scale(log(DIST)) + scale(log(LDIST)) + scale(log(AREA)) + 
##     fGRAZE + scale(ALT) + scale(YR.ISOL)
##  algorithm:    sampling
##  sample:       1800 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 56
##  predictors:   10
## 
## Estimates:
##                     mean   sd   10%   50%   90%
## (Intercept)        3.1    0.1  2.9   3.1   3.2 
## scale(log(DIST))   0.0    0.1 -0.1   0.0   0.1 
## scale(log(LDIST))  0.1    0.1  0.0   0.1   0.1 
## scale(log(AREA))   0.2    0.1  0.1   0.2   0.3 
## fGRAZE2            0.0    0.2 -0.2   0.0   0.2 
## fGRAZE3            0.0    0.1 -0.2   0.0   0.2 
## fGRAZE4            0.0    0.2 -0.2   0.0   0.2 
## fGRAZE5           -1.2    0.4 -1.8  -1.2  -0.7 
## scale(ALT)         0.0    0.1 -0.1   0.0   0.1 
## scale(YR.ISOL)     0.0    0.1 -0.1   0.0   0.1 
## sigma              6.6    0.7  5.7   6.6   7.5 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 19.2    1.2 17.7  19.2  20.8 
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##                   mcse Rhat n_eff
## (Intercept)       0.0  1.0  1621 
## scale(log(DIST))  0.0  1.0  1824 
## scale(log(LDIST)) 0.0  1.0  1813 
## scale(log(AREA))  0.0  1.0  1726 
## fGRAZE2           0.0  1.0  1872 
## fGRAZE3           0.0  1.0  1415 
## fGRAZE4           0.0  1.0  1629 
## fGRAZE5           0.0  1.0  1793 
## scale(ALT)        0.0  1.0  1792 
## scale(YR.ISOL)    0.0  1.0  1674 
## sigma             0.0  1.0  1568 
## mean_PPD          0.0  1.0  1819 
## log-posterior     0.1  1.0  1554 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Conclusions:

• in the Model Info, we are informed that the total MCMC posterior sample size is 1800 and that there were 56 raw observations.
• the estimated intercept (expected bird abundance when grazing intensity is equal to 1 and all of the continuous predictors are at their respective averages) is 3.06. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 21.4.
  
• the estimated slope associated with DIST (rate at which the abundance of birds changes per 1 unit change in log-transformed Distance to nearest patch holding other continuous predictors constant and grazing intensity at level 1), is -0.01 (mean) or -0.09 (median) with a standard deviation of 0. The 90% credibility intervals indicate that we are 90% confident that the slope is between -0.01 and -0.02 - e.g. there is no evidence of a significant trend.
  
• associated with each of the continuous predictors is a partial slope. Each partial slope is the rate of change between bird abundance and the associated predictor (on the log scale due to the link and based on 1 unit change in the predictor on the scale of the predictor). For example, for every one unit change in centred log patch Area, bird abundance is expected to increase by (log) 0.21.
• if we back transform this slope (by exponentiation), we get a partial slope for centred log Area of 1.23. This is interpreted as - for every 1 unit increase in (scaled log) Area, the bird abundance is expected to increase 1.23 fold. That is, there is a 23.08 percent increase per 1 unit increase in centred log Area.
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

tidyMCMC

tidyMCMC(loyn.rstanarm3$stanfit,
  estimate.method = "median", conf.int = TRUE,
  conf.method = "HPDinterval", rhat = TRUE, ess = TRUE
)

Conclusions:

• the estimated intercept (expected bird abundance when grazing intensity is equal to 1 and all of the continuous predictors are at their respective averages) is NA. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes NA.
  
• the estimated slope associated with DIST (rate at which the abundance of birds changes per 1 unit change in log-transformed Distance to nearest patch holding other continuous predictors constant and grazing intensity at level 1), is NA (mean) or -0.14 with a standard deviation of -0.02. The 90% credibility intervals indicate that we are 90% confident that the slope is between NA and 0.11 - e.g. there is no evidence of a significant trend.
  
• associated with each of the continuous predictors is a partial slope. Each partial slope is the rate of change between bird abundance and the associated predictor (on the log scale due to the link and based on 1 unit change in the predictor on the scale of the predictor). For example, for every one unit change in centred log patch Area, bird abundance is expected to increase by (log) NA.
• if we back transform this slope (by exponentiation), we get a partial slope for centred log Area of NA. This is interpreted as - for every 1 unit increase in (scaled log) Area, the bird abundance is expected to increase NA fold. That is, there is a NA percent increase per 1 unit increase in centred log Area.
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

gather_draws

loyn.rstanarm3 %>% get_variables()

##  [1] "(Intercept)"       "scale(log(DIST))"  "scale(log(LDIST))"
##  [4] "scale(log(AREA))"  "fGRAZE2"           "fGRAZE3"          
##  [7] "fGRAZE4"           "fGRAZE5"           "scale(ALT)"       
## [10] "scale(YR.ISOL)"    "sigma"             "accept_stat__"    
## [13] "stepsize__"        "treedepth__"       "n_leapfrog__"     
## [16] "divergent__"       "energy__"

loyn.draw <- loyn.rstanarm3 %>% gather_draws(`.Intercept.*|.*AREA.*|.*DIST.*|.*GRAZE.*|.*ALT.*|.*YR.*`, regex = TRUE)
loyn.draw

bayesplot

loyn.rstanarm3 %>% plot(plotfun = "mcmc_intervals")

tidy_draws

loyn.rstanarm3 %>% tidy_draws()

spread_draws

loyn.rstanarm3 %>% spread_draws(`.Intercept.*|.*DIST.*|.*AREA.*|.*GRAZE.*|.*ALT.*|.*YR.*`, regex = TRUE)

posterior_samples

loyn.rstanarm3 %>%
  posterior_samples() %>%
  as.tibble()

\(R^2\)

loyn.rstanarm3 %>%
  bayes_R2() %>%
  median_hdci()

brms

rstanarm captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

loyn.brm3 %>% summary()

##  Family: gaussian 
##   Links: mu = log; sigma = identity 
## Formula: ABUND ~ scale(log(DIST)) + scale(log(LDIST)) + scale(log(AREA)) + fGRAZE + scale(ALT) + scale(YR.ISOL) 
##    Data: loyn (Number of observations: 56) 
## Samples: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup samples = 2400
## 
## Population-Level Effects: 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept         3.04      0.11     2.82     3.25 1.00     2238     2208
## scalelogDIST     -0.01      0.06    -0.12     0.10 1.00     2150     2054
## scalelogLDIST     0.05      0.06    -0.07     0.17 1.00     2421     2277
## scalelogAREA      0.21      0.06     0.10     0.33 1.00     2321     2098
## fGRAZE2           0.03      0.15    -0.26     0.31 1.00     2237     2213
## fGRAZE3           0.05      0.14    -0.22     0.30 1.00     2159     2214
## fGRAZE4           0.01      0.16    -0.32     0.31 1.00     2296     2287
## fGRAZE5          -1.06      0.33    -1.79    -0.46 1.00     2648     2102
## scaleALT          0.00      0.05    -0.10     0.10 1.00     2279     2407
## scaleYR.ISOL      0.01      0.08    -0.13     0.16 1.00     2522     2497
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     6.30      0.63     5.20     7.74 1.00     2338     2364
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

• in the Model Info, we are informed that the total MCMC posterior sample size is 1800 and that there were 56 raw observations.
• the estimated intercept (expected bird abundance when grazing intensity is equal to 1 and all of the continuous predictors are at their respective averages) is 3.04. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 21.
  
• the estimated slope associated with DIST (rate at which the abundance of birds changes per 1 unit change in log-transformed Distance to nearest patch holding other continuous predictors constant and grazing intensity at level 1), is -0.01 (mean) or 0.1 (median) with a standard deviation of 0.06. The 90% credibility intervals indicate that we are 90% confident that the slope is between -0.01 and 1 - e.g. there is no evidence of a significant trend.
  
• associated with each of the continuous predictors is a partial slope. Each partial slope is the rate of change between bird abundance and the associated predictor (on the log scale due to the link and based on 1 unit change in the predictor on the scale of the predictor). For example, for every one unit change in centred log patch Area, bird abundance is expected to increase by (log) 0.21.
• if we back transform this slope (by exponentiation), we get a partial slope for centred log Area of 1.24. This is interpreted as - for every 1 unit increase in (scaled log) Area, the bird abundance is expected to increase 1.24 fold. That is, there is a 23.85 percent increase per 1 unit increase in centred log Area.
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

tidyMCMC

loyn.brm3$fit %>%
  tidyMCMC(
    estimate.method = "median",
    conf.int = TRUE, conf.method = "HPDinterval",
    rhat = TRUE, ess = TRUE
  )

Conclusions:

• the estimated intercept (expected bird abundance when grazing intensity is equal to 1 and all of the continuous predictors are at their respective averages) is NA. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes NA.
  
• the estimated slope associated with DIST (rate at which the abundance of birds changes per 1 unit change in log-transformed Distance to nearest patch holding other continuous predictors constant and grazing intensity at level 1), is NA (mean) or -0.13 with a standard deviation of -0.01. The 90% credibility intervals indicate that we are 90% confident that the slope is between NA and 0.1 - e.g. there is no evidence of a significant trend.
  
• associated with each of the continuous predictors is a partial slope. Each partial slope is the rate of change between bird abundance and the associated predictor (on the log scale due to the link and based on 1 unit change in the predictor on the scale of the predictor). For example, for every one unit change in centred log patch Area, bird abundance is expected to increase by (log) NA.
• if we back transform this slope (by exponentiation), we get a partial slope for centred log Area of NA. This is interpreted as - for every 1 unit increase in (scaled log) Area, the bird abundance is expected to increase NA fold. That is, there is a NA percent increase per 1 unit increase in centred log Area.
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

gather_draws

loyn.brm3 %>% get_variables()

##  [1] "b_Intercept"     "b_scalelogDIST"  "b_scalelogLDIST" "b_scalelogAREA" 
##  [5] "b_fGRAZE2"       "b_fGRAZE3"       "b_fGRAZE4"       "b_fGRAZE5"      
##  [9] "b_scaleALT"      "b_scaleYR.ISOL"  "sigma"           "prior_Intercept"
## [13] "prior_b"         "prior_sigma"     "lp__"            "accept_stat__"  
## [17] "stepsize__"      "treedepth__"     "n_leapfrog__"    "divergent__"    
## [21] "energy__"

loyn.draw <- loyn.brm3 %>% gather_draws(`^b_.*`, regex = TRUE)
loyn.draw

bayesplot

loyn.brm3 %>% mcmc_plot(type = "intervals")

tidy_draws

loyn.brm3 %>% tidy_draws()

spread_draws

loyn.brm3 %>% spread_draws(`^b_.*`, regex = TRUE)

posterior_samples

loyn.brm3 %>%
  posterior_samples() %>%
  as_tibble()

\(R^2\)

loyn.brm3 %>%
  bayes_R2(summary = FALSE) %>%
  median_hdci()

Further analyses

Model selection

In the frequentist instalment of this example, we explored three different approaches to pursuing parsimonious models. The first of these was to ‘dredge’ through all possible combinations of predictors and compare the resulting models via AICc. This approach has been widely criticised as it has the potential to yield models that are numerically ‘best’, yet ecologically meaningless.

The second approach was to average together the parameter estimates from many combinations of models so as to yield parameter estimates that are less dependent on the exact combination of predictors included in the model.

The final approach was to a priori propose a small (<10) number of possible candidate models, each of which focuses on a different aspect of the potential underlying ecological responses and see if any of those models are supported by the data. With this approach, the candidate models are likely to be sensible (as they are proposed on the basis of theory rather than purely numbers) and importantly, they are also interpretable. Furthermore, it also promotes the idea that multiple models might each offer important (and different) insights rather than there being a single ‘best’ model.

For this instalment, we will adopt the third approach. In the current context of birds in fragmented forests, we could propose the following models:

1. a model that focuses on the connectivity between patches (DIST and LDIST)
2. a model that focuses on the patch habitats (AREA, fGRAZE)
3. a model that focuses on the patch use and history (AREA, fGRAZE, YR.ISOL)
4. a model that focuses on climate (ALT)
5. a NULL model. The NULL model can be used as a reference point. Any model that has an AICc less (by at least 2 units) than (and thus better) than the NULL model must have some inferential support. Any model with AICc the same or higher than the NULL, clearly has no support and would not be considered a useful model

Note, these are all models that could be proposed prior to even collecting the data and all can be explained ecologically. By contrast, dredging by chance could suggest a model that has a combination of predictors that are very difficult to interpret and explain.

As the above models contain fewer predictors, there is now scope to include interactions.

Rather than use AIC, which is calculated once across the entire model, we can use other information criterion specifically designed for MCMC sampling. Of these, WAIC and loo are the most popular.

**Loo* (or leave-one-out cross validation based on the posterior likelihood), is actually an efficient approximation of leave-one-out cross validation. It does so using Pareto soothed importance sampling (PSIS).

• elpd_loo: expected log point-wise predictive density

• p_loo: effective number of parameters

• looic: LOO information criterion (-2 * elpd_loo)

• elpd_diff: returned by making pairwise comparisons between each model and the model with the largest ELPD (‘better model’) (which will be in the first row). We can put this in units of deviance by multiplying by -2

The reliability and approximate convergence rate of PSIS-based estimates can be inferred by examining the estimates of the Pareto distribution’s shape parameter (k):

• if k<0.5, the distribution of raw importance ratios has finite variance and central limit holds.
• if (0.5 <= k < 1), the distribution of raw importance ratios has infinite variance. Values under 0.7 are still useful, however values of k greater than 0.7 result in estimates that become unreliable. If there are k values between 0.7 and 1, moment matching should be performed.
• if k >= 1: bias in estimates is very large.

rstanarm

loyn.rstanarm4a <- update(loyn.rstanarm3, . ~ scale(log(DIST)) * scale(log(LDIST)),
  diagnostic_file = file.path(tempdir(), "dfa.csv")
)
loyn.rstanarm4b <- update(loyn.rstanarm3, . ~ scale(log(AREA)) * fGRAZE,
  diagnostic_file = file.path(tempdir(), "dfb.csv")
)
loyn.rstanarm4c <- update(loyn.rstanarm3, . ~ scale(log(AREA)) * fGRAZE * scale(YR.ISOL),
  diagnostic_file = file.path(tempdir(), "dfc.csv")
)
loyn.rstanarm4d <- update(loyn.rstanarm3, . ~ scale(ALT),
  diagnostic_file = file.path(tempdir(), "dfd.csv")
)
loyn.rstanarm4e <- update(loyn.rstanarm3, . ~ 1,
  diagnostic_file = file.path(tempdir(), "dfe.csv")
)
loo_compare(
  loo(loyn.rstanarm4a),
  loo(loyn.rstanarm4e)
)

##                 elpd_diff se_diff
## loyn.rstanarm4e  0.0       0.0   
## loyn.rstanarm4a -2.5       1.6

loo_compare(
  loo(loyn.rstanarm4b),
  loo(loyn.rstanarm4e)
)

##                 elpd_diff se_diff
## loyn.rstanarm4b   0.0       0.0  
## loyn.rstanarm4e -29.8       7.2

loo_compare(
  loo(loyn.rstanarm4c),
  loo(loyn.rstanarm4e)
)

##                 elpd_diff se_diff
## loyn.rstanarm4c   0.0       0.0  
## loyn.rstanarm4e -19.1       6.7

loo_compare(
  loo(loyn.rstanarm4d),
  loo(loyn.rstanarm4e)
)

##                 elpd_diff se_diff
## loyn.rstanarm4d  0.0       0.0   
## loyn.rstanarm4e -3.5       2.3

bayes_factor(
  bridge_sampler(loyn.rstanarm4a),
  bridge_sampler(loyn.rstanarm4e)
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

## Estimated Bayes factor in favor of x1 over x2: 0.00009

bayes_factor(
  bridge_sampler(loyn.rstanarm4b),
  bridge_sampler(loyn.rstanarm4e)
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

## Estimated Bayes factor in favor of x1 over x2: 803355.40456

bayes_factor(
  bridge_sampler(loyn.rstanarm4c),
  bridge_sampler(loyn.rstanarm4e)
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
## Iteration: 7
## Iteration: 8
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4

## Estimated Bayes factor in favor of x1 over x2: 0.01301

bayes_factor(
  bridge_sampler(loyn.rstanarm4d),
  bridge_sampler(loyn.rstanarm4e)
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

## Estimated Bayes factor in favor of x1 over x2: 1.99994

brm

loyn.brm4a <- update(loyn.brm3, . ~ scale(log(DIST)) * scale(log(LDIST)),
  save_pars = save_pars(all = TRUE), refresh = 0
)
loyn.brm4b <- update(loyn.brm3, . ~ scale(log(AREA)) * fGRAZE,
  save_pars = save_pars(all = TRUE), refresh = 0
)
loyn.brm4c <- update(loyn.brm3, . ~ scale(log(AREA)) * fGRAZE * scale(YR.ISOL),
  save_pars = save_pars(all = TRUE), refresh = 0
)
loyn.brm4d <- update(loyn.brm3, . ~ scale(ALT),
  save_pars = save_pars(all = TRUE), refresh = 0
)
loyn.brm4e <- update(loyn.brm3, . ~ 1,
  save_pars = save_pars(all = TRUE), refresh = 0
)
waic(loyn.brm4a)

## 
## Computed from 2400 by 56 log-likelihood matrix
## 
##           Estimate  SE
## elpd_waic   -216.2 4.0
## p_waic         5.0 0.9
## waic         432.5 7.9
## 
## 3 (5.4%) p_waic estimates greater than 0.4. We recommend trying loo instead.

loo(loyn.brm4a)

## 
## Computed from 2400 by 56 log-likelihood matrix
## 
##          Estimate  SE
## elpd_loo   -216.4 4.0
## p_loo         5.2 0.9
## looic       432.8 8.0
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     54    96.4%   874       
##  (0.5, 0.7]   (ok)        1     1.8%   850       
##    (0.7, 1]   (bad)       1     1.8%   270       
##    (1, Inf)   (very bad)  0     0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.

loo_compare(
  loo(loyn.brm4a),
  loo(loyn.brm4e)
)

##            elpd_diff se_diff
## loyn.brm4e  0.0       0.0   
## loyn.brm4a -2.5       1.7

# -2 * -2.5
loo_compare(
  loo(loyn.brm4b),
  loo(loyn.brm4e)
)

##            elpd_diff se_diff
## loyn.brm4b   0.0       0.0  
## loyn.brm4e -30.9       7.8

loo_compare(
  loo(loyn.brm4b, moment_match = TRUE),
  loo(loyn.brm4e)
)

##            elpd_diff se_diff
## loyn.brm4b   0.0       0.0  
## loyn.brm4e -30.8       7.8

loo_compare(
  loo(loyn.brm4c, moment_match = TRUE),
  loo(loyn.brm4e)
)

##            elpd_diff se_diff
## loyn.brm4c   0.0       0.0  
## loyn.brm4e -21.1       7.5

loo_compare(
  loo(loyn.brm4d),
  loo(loyn.brm4e)
)

##            elpd_diff se_diff
## loyn.brm4d  0.0       0.0   
## loyn.brm4e -3.7       2.6

An alternative is to compute Bayes factors based on bridge sampling. Note, this process usually requires far greater number of posterior samples (10x) in order to be stable. It is also advisable to run this multiple times to ensure stability.

It calculates the marginal likelihood of one model in favour of another. The larger the value, the more evidence there is of one model over the other.

bayes_factor(
  loyn.brm4a,
  loyn.brm4e
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

## Estimated Bayes factor in favor of loyn.brm4a over loyn.brm4e: 0.00129

# OR
bayes_factor(
  loyn.brm4e,
  loyn.brm4a
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4

## Estimated Bayes factor in favor of loyn.brm4e over loyn.brm4a: 769.48559

bayes_factor(
  loyn.brm4b,
  loyn.brm4e
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
## Iteration: 7
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

## Estimated Bayes factor in favor of loyn.brm4b over loyn.brm4e: 23178980133.56327

# OR
bayes_factor(
  loyn.brm4e,
  loyn.brm4b
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
## Iteration: 7

## Estimated Bayes factor in favor of loyn.brm4e over loyn.brm4b: 0.00000

bayes_factor(
  loyn.brm4c,
  loyn.brm4e
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

## Estimated Bayes factor in favor of loyn.brm4c over loyn.brm4e: 146839.98926

# OR
bayes_factor(
  loyn.brm4e,
  loyn.brm4c
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

## Estimated Bayes factor in favor of loyn.brm4e over loyn.brm4c: 0.00001

bayes_factor(
  loyn.brm4d,
  loyn.brm4e
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

## Estimated Bayes factor in favor of loyn.brm4d over loyn.brm4e: 8.17588

# OR
bayes_factor(
  loyn.brm4e,
  loyn.brm4d
)

## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

## Estimated Bayes factor in favor of loyn.brm4e over loyn.brm4d: 0.12253

Contrasts

brms

loyn.list <- with(loyn, list(AREA = c(min(AREA), mean(AREA), max(AREA))))

newdata <- loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  pairs() %>%
  as.data.frame()
head(newdata)

loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list) %>%
  regrid() %>%
  pairs()

## AREA =    0.1:
##  contrast estimate lower.HPD upper.HPD
##  1 - 2       9.294     0.694     19.11
##  1 - 3      10.323     0.894     19.59
##  1 - 4      16.191     7.431     25.33
##  1 - 5      17.065     9.666     25.87
##  2 - 3       0.966    -6.671      9.07
##  2 - 4       6.870    -0.534     14.44
##  2 - 5       7.736     1.068     14.53
##  3 - 4       5.915    -1.843     12.67
##  3 - 5       6.682     1.288     13.77
##  4 - 5       0.722    -3.801      7.83
## 
## AREA =   69.3:
##  contrast estimate lower.HPD upper.HPD
##  1 - 2      -3.183   -10.440      4.18
##  1 - 3      -5.365   -13.947      4.42
##  1 - 4     -24.992   -71.614      1.07
##  1 - 5       8.318   -11.813     23.73
##  2 - 3      -2.187   -13.059      9.20
##  2 - 4     -21.617   -66.587      7.96
##  2 - 5      11.241   -10.116     29.02
##  3 - 4     -19.818   -64.696     12.25
##  3 - 5      13.350    -9.618     30.33
##  4 - 5      32.628    -9.247     78.21
## 
## AREA = 1771.0:
##  contrast estimate lower.HPD upper.HPD
##  1 - 2     -21.710   -57.156      5.77
##  1 - 3     -30.583   -75.732      9.54
##  1 - 4    -203.249 -1151.671     10.59
##  1 - 5     -31.248  -308.030     36.73
##  2 - 3      -9.001   -62.776     46.04
##  2 - 4    -178.997 -1144.401     46.68
##  2 - 5     -10.066  -286.823    102.55
##  3 - 4    -171.191 -1140.089     68.90
##  3 - 5      -0.977  -282.146    125.12
##  4 - 5     155.627  -356.586   1343.55
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

newdata <- loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  pairs() %>%
  gather_emmeans_draws()

newdata %>%
  median_hdci() %>%
  ggplot() +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_pointrange(aes(y = .value, ymin = .lower, ymax = .upper, x = contrast)) +
  facet_wrap(~AREA) +
  coord_flip()

emmeans(loyn.brm4b, ~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  gather_emmeans_draws()

newdata.p <- newdata %>% summarise(P = sum(.value > 1) / n())
g <- newdata %>%
  ggplot() +
  geom_vline(xintercept = 1, linetype = "dashed") +
  stat_slab(aes(
    x = .value, y = contrast,
    fill = stat(ggdist::cut_cdf_qi(cdf,
      .width = c(0.5, 0.8, 0.95),
      labels = scales::percent_format()
    ))
  ), color = "black") +
  scale_fill_brewer("Interval", direction = -1, na.translate = FALSE) +
  facet_grid(~ round(AREA, 1))

g + geom_text(data = newdata.p, aes(y = contrast, x = 1, label = round(P, 3)))

g + geom_text(data = newdata.p, aes(y = contrast, x = 1, label = paste("P = ", round(P, 3))), hjust = -0.2, position = position_nudge(y = 0.5))

Summary figure

rstanarm

From the full model

loyn.list <- with(loyn, list(AREA = modelr::seq_range(AREA, n = 100)))

newdata <- emmeans(loyn.brm3, ~ AREA | fGRAZE, at = loyn.list, type = "response") %>%
  as.data.frame()
head(newdata)

ggplot(newdata, aes(y = response, x = AREA)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = fGRAZE), alpha = 0.3) +
  geom_line(aes(color = fGRAZE)) +
  theme_bw() +
  scale_x_log10() +
  scale_y_log10()

spaghetti <- emmeans(loyn.rstanarm3, ~ AREA | fGRAZE, at = loyn.list, type = "response") %>%
  gather_emmeans_draws() %>%
  mutate(Fit = exp(.value))
wch <- sample(1:max(spaghetti$.draw), 100, replace = FALSE)
spaghetti <- spaghetti %>%
  filter(.draw %in% wch)
ggplot(newdata) +
  geom_line(data = spaghetti, aes(
    y = Fit, x = AREA, color = fGRAZE,
    group = interaction(fGRAZE, .draw)
  ), alpha = 0.05) +
  geom_line(aes(y = response, x = AREA, color = fGRAZE)) +
  theme_bw() +
  scale_x_log10() +
  scale_y_log10()

From the “habitat” model

loyn.list <- with(loyn, list(AREA = seq(min(AREA), max(AREA), len = 100)))

newdata <- emmeans(loyn.rstanarm4b, ~ AREA | fGRAZE, at = loyn.list, type = "response") %>%
  as.data.frame()
head(newdata)

ggplot(newdata, aes(y = response, x = AREA)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = fGRAZE), alpha = 0.3) +
  geom_line(aes(color = fGRAZE)) +
  theme_bw() +
  scale_x_log10() +
  scale_y_log10()

spaghetti <- emmeans(loyn.rstanarm4b, ~ AREA | fGRAZE, at = loyn.list, type = "response") %>%
  gather_emmeans_draws() %>%
  mutate(Fit = exp(.value))
wch <- sample(1:max(spaghetti$.draw), 100, replace = FALSE)
spaghetti <- spaghetti %>% filter(.draw %in% wch)
ggplot(newdata) +
  geom_line(data = spaghetti, aes(
    y = Fit, x = AREA, color = fGRAZE,
    group = interaction(fGRAZE, .draw)
  ), alpha = 0.05) +
  geom_line(aes(y = response, x = AREA, color = fGRAZE)) +
  theme_bw() +
  scale_x_log10() +
  scale_y_log10()

brms

From the full model

loyn.list <- with(loyn, list(AREA = modelr::seq_range(AREA, n = 100)))

newdata <- emmeans(loyn.brm3, ~ AREA | fGRAZE, at = loyn.list, type = "response") %>%
  as.data.frame()
head(newdata)

ggplot(newdata, aes(y = response, x = AREA)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = fGRAZE), alpha = 0.3) +
  geom_line(aes(color = fGRAZE)) +
  theme_bw() +
  scale_x_log10() +
  scale_y_log10()

spaghetti <- emmeans(loyn.brm3, ~ AREA | fGRAZE, at = loyn.list, type = "response") %>%
  gather_emmeans_draws() %>%
  mutate(Fit = exp(.value))
wch <- sample(1:max(spaghetti$.draw), 100, replace = FALSE)
spaghetti <- spaghetti %>% filter(.draw %in% wch)
ggplot(newdata) +
  geom_line(data = spaghetti, aes(
    y = Fit, x = AREA, color = fGRAZE,
    group = interaction(fGRAZE, .draw)
  ), alpha = 0.1) +
  geom_line(aes(y = response, x = AREA, color = fGRAZE)) +
  theme_bw() +
  scale_x_log10() +
  scale_y_log10()

From the “habitat” model

Lets explore the relationship between bird abundance and patch area for each grazing intensity separately.

loyn.list <- with(loyn, list(AREA = modelr::seq_range(AREA, n = 100)))

newdata <- emmeans(loyn.brm4b, ~ AREA | fGRAZE, at = loyn.list, type = "response") %>%
  as.data.frame()
head(newdata)

ggplot(newdata, aes(y = response, x = AREA)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = fGRAZE), alpha = 0.3) +
  geom_line(aes(color = fGRAZE)) +
  theme_bw() +
  scale_x_log10() +
  scale_y_log10()

spaghetti <- emmeans(loyn.brm4b, ~ AREA | fGRAZE, at = loyn.list, type = "response") %>%
  gather_emmeans_draws() %>%
  mutate(Fit = exp(.value))
wch <- sample(1:max(spaghetti$.draw), 100, replace = FALSE)
spaghetti <- spaghetti %>% filter(.draw %in% wch)
ggplot(newdata) +
  geom_line(data = spaghetti, aes(
    y = Fit, x = AREA, color = fGRAZE,
    group = interaction(fGRAZE, .draw)
  ), alpha = 0.1) +
  geom_line(aes(y = response, x = AREA, color = fGRAZE)) +
  theme_bw() +
  scale_x_log10() +
  scale_y_log10()

References

Loyn, R. H. 1987. “Nature Conservation: The Role of Remnants of Native Vegetation.” In, edited by D. A. Saunders, G. W. Arnold, A. A. Burbridge, and A. J. M. Hopkins. Chipping Norton, NSW: Surrey Beatty & Sons.